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Abstract 

Conditional Asian options are recent market innovations, which offer cheaper and long-dated 
alternatives to regular Asian options. In contrast with payoffs from regular Asian options which 
are based on average asset prices, the payoffs from conditional Asian options are determined only 
by average prices above certain threshold. Due to the limited inclusion of prices, conditional 
Asian options further reduce the volatility in the payoffs than their regular counterparts and have 
been promoted in the market as viable hedging and risk management instruments for equity- 
linked life insurance products. There has been no previous academic literature on this subject 
and practitioners have only been known to price these products by simulations. We propose 
the first analytical approach to computing prices and deltas of conditional Asian options in 
comparison with regular Asian options. In the numerical examples, we put to the test some 
cost-benefit claims by practitioners. As a by-product, the work also presents some distributional 
properties of the occupation time and the time-integral of geometric Brownian motion during 
the occupation time. 

Key Words. Option pricing; hedging; conditional Asian option; Asian option; Laplace 
transform inversion; asymptotic expansion; integral of geometric Brownian motion; occupation 
time; Lommel functions. 


1 Introduction 


Asian options (average options) are widely used in commodity and stock markets as cheaper al¬ 
ternatives to European and American options for hedging and risk management. A common use 
of Asian options is to transfer the risk of volatility in asset prices to capital markets. Consider 
a continuous-time stochastic model for asset prices, denoted by {Xf,t > 0} with Xq = x. The 
payoff from a continuous European-style Asian put option with fixed strike K and maturity T is 
determined by 



( 1 ) 


where (x)_|_ = max{x,0}. A review of historical references and examples of Asian options can be 
found in Linetsky, 2004 . The pricing of Asian options led to a flourishing of computational tech¬ 
niques in the mathematical finance literature. Monte Carlo simulations were among the first pricing 
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methods, proposed by |Kemna and Vorst, 1992| , |Boyle et ah, 1997| . Bounds and approximations 
based on conditional expectations can be seen in |Curran, 1994| , |Nielsen and Sandmann,~2 002|, 

[Lord, 2006]. Numerical PDE methods were extensively studied in a series of papers including 
|Rogers and Shi, 1995| , |Alziary et ah, 1997] , |Vecer, 2001| , |Zvan et al., 2000| . Numerical inver¬ 
sion of Laplace/Fourier transforms were developed in |Geman and Yor, 1993| , |Carverhill and Clewlow,~1 9901, 
|Fu et al., 1998| , |Cai and Kou, 2012] , |Carr and Schroder, 2003 . Latest work by |Chen et ah, 2012] 
also extends Hilbert transform methods to Asian option pricing. A variety of orthogonal poly¬ 
nomial/eigenfunction expansions were used in |Dufresne, 2000| , |Linetsky, 2004| , |Ju, 200 21, and 
|Zhang and Oosterlee, 2013 


Asymptotics expansions of Asian option prices were developed in 


et al., 1997| ; |Cai et ah, 2014| , |Gerhold, 2011 , etc. Given the sheer volume of research 


[Carmona 

papers on this subject matter, the list of aforementioned papers is by no means a comprehensive 
review of the literature. 

In recent years, BNP Paribas introduced a variation of the Asian option, termed conditional 
Asian option, under which the average of asset prices is only based on prices that are above a 
certain threshold (See an industrial presentation by |Segau d, 20lT]). Formulated mathematically, 
the payoff from the continuous conditional Asian option is given by 


K — 


f 5 

Jo 


Xtl{x t >bj dt 


f 

Jo 


J{Xt>b} d t 


(2) 


where the threshold level is denoted by b > 0 and I a is an indicator function which equals 1 if A 
is true and 0 otherwise. Figure |T] gives a sample path of asset prices (black solid line) where the 
average-to-date (red dash line) only includes prices above a threshold (blue dotted line). Observe 
that the sample path of the average-to-date becomes flat whenever asset prices plunge below the 
threshold. Parameters used in this numerical example are provided in Section [9j In the BNP 



Figure 1: Example of average-to-date above a threshold. 

Paribas presentation, the threshold b is called observation barrier and is 50% of the initial price xq. 
They gave an illustrative example of a five-year at-the-money conditional Asian option (K = Xq), 
which is the discrete version of (O with monthly observations. 
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There are several reasons why conditional Asian puts may be considered effective alternatives 
to regular Asian puts in the financial market. First, conditional Asian puts are cheaper than Asian 
puts with the same term and strike. Observe that 
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X t dt. 


Thus the payoff from a conditional Asian option Q is surely no greater than that from a reg¬ 
ular Asian option (HD- According to BNP Paribas, a five-year at-the-money conditional Asian 
option was shown to cost 40% less in comparison with a regular Asian option, yet achieved 75% 
of the initial delta for a regular one. Second, most Asian options available in the market are rel¬ 
atively short termed (typically less than one year). In fact, there are very few hedging products 
that can match the duration of insurance liabilities due to their long-term nature. BNP Paribas 
promoted five-year conditional Asian options as hedging solutions for variable annuities, which 
are investment-combined insurance products with very long term (typically more than ten years). 
Third, cash flow structures of variable annuity guaranteed benefits resemble that of Asian options. 
This observation was exploited in |Feng and Volkmer, 2012} Feng and Volkmer, 2014, Feng, 2014 


for the computation of risk measures for guaranteed benefits. It follows from the same logic that 
long-term Asian-like options provide natural hedge against financial risks embedded in variable 
annuity products. Therefore, conditional Asian options are touted as a cost effective tool for asset 
and liability management. 

In this paper, we intend to develop an analytical method to price the conditional Asian option 
and compute its delta. We consider the Black-Scholes model under which the underlying asset 
follows the geometric Brownian motion under the risk-neutral measure 


dX t = rX t dt + aX t dB t , 


and under which the no-arbitrage price of a conditional Asian put option can be written as 


AP h :=E 


-rT 


K — 


lXJ{Xt>b} dt 

fo I {X t >bj dt , 


where E denotes the expectation under the risk-neutral probability measure. The regular Asian 
put is clearly a special case of conditional Asian put with 6 = 0 and hence its price is given by APo- 
It turns out that the price of conditional Asian put is best computed as the price of regular Asian, 
which can be done through inversion of a closed-form Laplace transform, minus some price spread, 
which needs to be worked out by inversion of a Laplace transform with an integral representation. 


AP i) = APq — Price Spread. 


The price spread is clearly a result of cost saving from filtering out undesirable low prices in the 
average price. A practical advantage of such a split is that one can compute the price spread directly 
without having to carry out simulations for pricing conditional Asian options. Practitioners may 
use it to extrapolate conditional Asian prices by combining market prices of Asian options with 
estimated price spread. The deterministic algorithm developed later for the computation of price 
spread can be easily used with implied volatility, making the conditional Asian option pricing 
consistent with the standard practice of pricing option with less liquidity and thin trading. 
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2 Asian options revisited 


As it shall become clear later, the regular Asian option plays an important role in the pricing of 
conditional Asian options. We shall present some results on Asian put option prices, which can be 
shown to be consistent to known results on Asian call option prices in the literature. 

Consider the process 


Y t ■= ( X s ds, t> 0. 

Jo 

Let p(x, t, y) be the transition density of Y with Xq = x and set 


rv 

P{x,t,y)= p(x,t,v)dv, Q(x, t, y) 
Jo 



P{x, t, v ) dv. 


Then it follows immediately that 


(3) 

(4) 


AP 0 = -e~ rl Q(x,T,TK), 


(5) 


In what follows, we derive an explicit solution to the Laplace transform of Q(x,t,y) with respect 
to t only. This allows computation of Q by a one-dimensional numerical Laplace inversion. 

In |Feng and Volkmer, 2014] we considered the special case 

X f = exp(2z/t + 2B t ), 2 )t = / X s d s. 

Jo 

Setting u = Aw /a 2 and using the scaling property B ct ~ y/cB t , we get 


Yt = / X u du~x 


where the constant v is given by 


± 

v 2 Jo 


Ax , 


exp(2 uw + 2 B w ) dw = — 2 } CT 2 t / 4 , 


(T 


2r 

u=--l. 

cr A 

We denote the transition density of 2) by po and define P $, Q$ in the same way in which P , Q are 
defined by p in dH). It follows immediately that 


P(x, t, y) = P 0 


a 2 t a 2 y 

"T’laT 


( 6 ) 


cr 2 fcr 2 ta 2 y 


Ax „ f a 2 t a 2 y 


Q{x, t, y) = —jQo —t~ i —7— 


(T 


A ’ Ax 


( 7 ) 


Therefore, the remaining task is to find a way of computing Qq. 

We shall denote the Laplace transform f(s) := e~ st f(t) df. It is known that 


Po(s,y) := f e st p 0 (t,y)dt.= ^ K A(y) 5 

Jo r(l + 2 rj) 


( 8 ) 
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where 


f K (y) = y K exp ( M K>?7 (J- j , 


with M k ^ denoting a Whittaker function, and 


1 — v 1 /- 

k =—-—, r/ = -y2s + 


We use integration by parts to show that 




It follows that 


Pq (s,y) = 


r (j]-k + \) 


2(77 + k — A) (77 - k + A) (77 + k — A)T(1 + 2rj) 


2~ K L-i(y) 


(9) 


and 


Qo(s,y) = 


2(77 + K - \){rj - K + \) 4(77 + K - \){r) + k - |)(77 - K + |)(?7 - K + 5) 


+ 


r(?7- « + i) 


F(1 + 277 ) (77 + re-i )(77 + /e-§) 


fn- 2 (y)- 


( 10 ) 

Although the expression for Qo was n °t explicitly stated in the literature, similar expressions are 
known for Asian call options. |Gen ian and Yor, 1993, (3.10)] developed a single integral expression 
for a similar Laplace transform of Asian call price with respect to the time parameter, which was 
later improved by Donati-Martin et ah, 2001 in terms of Kummer’s function of the first kind and 
similarly by [Shaw, 2003, Section 3]. It can be shown through the put-call parity that the result 
here is consistent with that of Asian call price in Donati-Martin et ah, 2001, (3.8)]. In the first 
example of Section [9] we shall verify numerically the equivalence of the formula in this section and 
those known in the literature. 


3 Joint Laplace transform 


The computation of the conditional Asian option price requires the study of the joint distribution 
of the pair ( Ut,Vt ) where 


L 


U t = I{x T >b} dr, V t — X T I{x T >b}dT, 


f 

Jo 


b> 0. 


The former is known as an occupation time whose transition density is studied in |Pcchtl, 1999 
The latter is known in the literature only in the case of b = 0, which is also called Yor’s process. 
Let us consider the triple joint Laplace transform 


F(b, x, s, a, /3) 


^t x [exp{-aU Ts -PV Ts }], 


( 11 ) 
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where T s is an exponential random variable with mean 1/s independent of X. By the Feynman-Kac 
theorem, F(x ) = F(b, x, s, a, f3) satisfies the differential equation 

\a 2 x 2 F"(x) + rxF'(x) — [(a + /3x)I{ x> b} + s] F(x) = —1, x > 0 (12) 

with boundary conditions 

lim F(x) = 

2-S-0+ s 

lim F(x) = 0. (13) 

x —^OO 

As usual, a solution F(x) of (fT2l) is understood to be differentiable at x = b with a jump in the 
second derivative there. 


3.1 The case 6 = 0 

We first consider the simpler case where 6 = 0. In this case, Ut = t and 

F(0, x,s,a,(3) = -E'/expj -aT s - /3Vr s }\. 
s 

Then F(x) = F(0, x, s,a, /3) satisfies the ODE 

\< j 2 x 2 F” ( x ) + rxF’{x) — (a + s + /3x)F(x) = —1, x > 0, 
subject to the boundary conditions (fl3l) and 

lim F(x) = -. 

ai-S"0+ s + a 


The homogeneous equation 

\o 2 x 2 F" (x) + rxF'(x) — (a + s + (3x)F(x) = 0 
has the fundamental system of solutions 

VWxj , 

F 2 {x)=x~^ l+l1 ^ 2 K\ \f2f3x 

where I\,K\ denote modified Bessel functions and 

2r 

d 9 2 , 

a z 

A= \f{n + l) 2 + 8 a~ 2 (s + a). 



F\(x) = x ( 1+ ^)/ 2 / a 


(14) 


(15) 


(16) 


(17) 


When a,(3,s > 0, the solution F\ (x) is positive and increasing on (0,oo) with linx E _>.o+ Fi (%) = 0, 
lim^^oo F\ (x) = oo while F 2 (x) is positive and decreasing on (0,oo) with lim x ^o+ F 2 {x) = oo, 
lim^^oo F 2 (x ) = 0. The Wronskian is 

Fi(x)F'(x) - F[(x)F 2 (x) = -\x~ 2 -». 
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From the variation of parameters formula, we obtain 


a poo 4 r x 

F(x) = —zF 1 (x) z^F 2 (z)dz + —F 2 (x) z^Fi(z)dz. (18) 

J x Jo 

Note that the integrals in (1181) are well-defined. 

Certainly F(x) defined by (fT8j) satisfies (fTTl) . To verify (fl3l) we use the known behavior of 
K\(z) as z —>• 0: 


h(z)' 


(§) ; 


r(A + i)’ 


x»W~ir W (f) 
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Then, by using L’Hospital’s rule, we get 

4 1 
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1 


i I o/( ;1 ' ) a 2 A(A — /x — 1) ' a 2 A(A + // + 1) 
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s + a 


(19) 

( 20 ) 


In a similar way, we can show that also (USD is true. Note that since lim^oo F\ (x) = oo and 
lim x ^o+ F 2 (x) = oo, there is only one solution of (fT4]) . (fl3l) . (fl5l) . so F{x) must be the function we 
are looking for. 

In the Appendix we use the theory of Lommel’s functions to show that F in (|18D can be written 
in the form 


1 


F( 0, x, s,a,/3) = -iF 2 (1; Ufi - A + 3), U/i + A + 3); 2cr 2 fix) 


s + a 


+ ~ ^ + l))r(3(M + ^ + 1)) Fi(x ), (21) 


where \F 2 is a generalized hypergeometric function defined for all z € C by the convergent series 


1 -^ 2 (&i! bi, b 2 ] z) — y 
k =0 


(oi )_k z k 

(h )k(b 2 ) k k\ ’ 


{biM ^ 0,-1,—2,-■ ■) 


with the Pochhammer’s symbol ( a) n = a(a + 1) ■ ■ ■ (a + n — 1) if n = 1, 2, • ■ ■ and (a)o = 1. 

The advantage of (1211) over (fl8l) is that it does not involve an integration. The disadvantage is 
that (|2l1) breaks down when A — ^ = 2^+1,^ = 1,2,... The condition A — /r = 2p + 1 is equivalent 
to 



2(s + a) 


+ p{p - 1) = 0. 


However, we will use m only for purely imaginary values of a,/3 and then the exceptional cases 
do not occur. Another problem with (1211) is that we might have huge cancellations in the addition 
of the two summands. Using arbitrary precision software the problem can be solved by increasing 
the number of digits although this is not an elegant method. Later we will make sure that we use 
m only in cases when there are no large cancellations. 
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3.2 The case b > 0 

If x > b then there is a constant A such that 


F(x) = AF 2 (x) + Y(x), 

where F2 is defined as in the previous subsection and Y is given by (1181) or m, i-e. 

Y(x) = F( 0, x, s, a, (3). 

Note that we could choose Y (x) differently, but it must be a solution of (|12l) which converges to 0 
as x —>• 00 . If 0 < x < b, then there is a constant B such that 

F(x) = Bx p + -, 
s 

where 


P— + 1) + 5 ^ 0 , 

A 0 = y/(n + l) 2 + 8a~ 2 s. 

By matching the values of F and F' at x = b we obtain for x > b 


F(b, x, s, a, /3) = d>( 6 , x, s, a, f3) + Y(x), 


where 


<h( 6 , x, s, a, /3) 


p(\-Y(b))+bY'(b) 

P F 2 (b)-bF'(b) 


F 2 {x). 


Observe that as b —>• 0, $(6, x, s, a, f3) approaches zero due to (fT3|) . 


( 22 ) 


(23) 


4 Moments of (Ut, Vt) 

We observe that when a = irz and (3 = —it, the mapping r —>• sF(b, x, s, a, (3) becomes the 
characteristic function of W := Vr a — zUt s , i-e. 

sF(b,x, s, irz,— it) = E x [exp{irIT}]. 

In this section we will determine the coefficients Ek in the expansion 


sF(b, x, s, iTZ, —it) = Eq + E\t + • • • + E p t p + o(r p ) as r —>• 0. (24) 

Then the moments of W are given by 

E x [W k ] = i~ k k\E k \ 

see |Lukacs, 1970[ Theorem 2.3.3]. Not only does (1M1) lead to a way of computing moments of 
( Ut , Vt), it also shows the rate of convergence near r = 0 of an important integral of F to be seen 
in the computation of conditional Asian option prices. 

Consider the case where x > b. Let p € N. We assume that p > p which is equivalent to 

2s 2 r 

— >p{p- 1) +p —• 


(7 










For example, if p = 1 then we require that s > r. We set a = irz, (3 = —It in m and use (USD. 
Since p > p we see that the second summand on the right-hand side of ( 1211 ) is o(t p ) as r — )• 0. The 
first term on the right-hand side of (1211) is analytic at r = 0. Therefore, we can write 

sY ( x , s, irz, —it) = Aq + A\t + • • • + A p t p + o(r p ), 


where Aq = 1 and 


A\ = —i 


x z 

-1— 

r — s s 


In a similar way we get (' = d/dx) 

sY'(x, s, irz, —it) = B\t + B 2 T 2 + • • • + B p r p + o(r p ), 

where 


B\ = — i- 


1 


r — s 


It follows that 

where 


p(l - sY(b)) + bsY\b) = Cir + • • • + C p t p + o{t p ), 


Ci = ip 


b z\ b 


r — s s J r — s 
Now consider the function (recall a = irz , /3 = —ir ) 

= F 2 (x) 

pF 2 (b)-bF>(bY 

We express uj in terms of the Bessel function I\ and I-\ by using (fTP) , 

7 r 


(25) 


K v (u) — . AI- V (u) I v (u)) 

Z sm(7rx) 


and 




u 


2k 


k =0 


4: k k\T(v + k + 1) 


Upon substitution of K\ in (fT71) . we observe that the term involving I\ is 0(/3 A / 2 ) as t —» 0. Thus, 

F 2 (x) = / x-^/^f-v^") + O (> /2 ) . 

2sm(A7r) \cr J \ J 

We treat the term pF 2 (b ) — bF^b) similarly: 
pF 2 (» - bFl(b) = hzAb-«+^ Kx v^j) + Mb-^K M (Zp^b) 

* (Ijw) - ' + o b 

\cr J 2sin(A7r) a \a J \ 


A/2 


2sin(A7r) 2 


Multiplying the numerator and the denominator of the fraction in (1251) by ^ sin(7rA) T(—A) 

gives 


“M = (|) 


-(l+/i+A)/2 ^ (Jj./Tc)*' 


k=0 


4 fc fc!(-A) fc+1 


E 

Uc=0 


^o_+A _ ^ (^-/3 & ) 


k 1 - 1 


4*fc!(-A) fc+1 


+ 0(r p ), (26) 
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where we used that Ao + p + 1 > 0 implies \q > p > p. Expanding the right-hand side in powers 
of r yields 

ui(t) = D 0 + Dpr -I-+ D p t p + o(t p ), 


where 


A) = T- ( w 

An 


When we put everything together, we obtain 


— 5OH-A0+I) 

with Eq = 1 and 
x\ -^hH-Ao+i) 


— i 


r — s 


Similarly, one can also work out the expansion when 0 < x < b and 


Ei = —i 


b z 

-1— 

r — s s 


/r + Aq + 1 


+ 


1 fx\P 


r — s I Aq V b / 



Therefore, 


lim — i sF(b,x,s,iTZ,—iT) = — 
t—> o t i 


By identifying the coefficients of z’s, we obtain that when x > b, 

b{p — 1) fx\ -lO+Ao+i) x 


E(V Ts ) = 


('r - s)A 0 \b 


r — s 


S AqS 


(27) 


and when 0 < x < b, 


Efe) 

K(U Ts ) 


b(p + A 0 - 3) /x\p 
2(s — r)Ao 
/r + Ao + 1 / x\p 
2Aq s \b) 


The above method can be used to evaluate the coefficients Ek for k > 2 with any computer 
algebra system, such as Maple and Mathematica, and hence higher moments and cross-moments 
of ( Ut s , Vt s )- We have developed a computer routine to compute these coefficients, which shall be 
made available on the authors’ webpages. 


5 Conditional Asian option 

The key element of pricing conditional Asian option is to find the distribution of the ratio 

V t 

Z t := —, t > 0. 

U t 

We denote its cumulative distribution function by 

G{b , x , z, t ) = P (Z t < z) 

Since Vt/Ut > 6, we must have 

G(b, x, z,t) = 0 for z <b. 


(28) 


(29) 
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To price a conditional Asian option with the payoff ([2]), we need to calculate 


AP b = e~ rT E((K - Z T )+) 



x, z, T) dz, 


(30) 


which can be shown through integration by parts. 

Applying Gurland’s formula 
Ut s /Vt s , we obtain for s > 0, 


Gurland, 1948, Theorem 1] to the ratio of random variables 


s 



e st G(b, x, z,t)dt 


11 1 

-P.V. / —sF(b,x,s,iTZ,—iT)dT. 

2 2iri J r 


By the reflection principle, we must have F(b. x, s, a, f$) = F(b, x, s, a, j3) and thus 



e st G(b,x, z,t) dt 


1 1 f°° 1 

-/ —2sF(b,x,s,iTZ,—iT) dr. 

2s ttJo t 


(31) 


As it is difficult to develop a closed-form solution to G itself, we want to carry out the compu¬ 
tation of the conditional Asian option price in (|30l) in three steps. First, we compute the integral 
on the right-hand side of (13111 numerically. In the previous section, we discussed the behavior of 
the integrand as t — >• 0 and we shall analyze its behavior as r —>• +oo. Then, we will use the 
Gaver-Stehfest algorithm for numerical inversion of Laplace transform to obtain G. Finally, we 
perform the numerical integration in (13011 . It is interesting to note that the three operations 

• Gurland integral (13T1) . 


• inverse Laplace transform s t, 

• integration of distribution function (13011 
could be carried out in any order. 

Remark 1. It is worthwhile to mention that the computation of the regular Asian option price using 
[Yor, 1992]’s integral representation of the joint density of the geometric Brownian motion and its 
time integral also requires the evaluation of a triple integral. Such computation was demonstrated in 
jls hiyama, 2 005]. Even though the payoff of conditional Asian option is considerably more evolved 
than its regular counterpart, the complexity of computation appears to be comparable. 

Remark 2. In the case ofb = 0 (i.e. regular Asian option), the computation of the Gurland integral 
(Em was avoided in Section [H because we exploited the fact that the two-dimensional random vector 
( Ut,Vt ) degenerates to the vector ( t,Y t ) driven by a single random variable and hence 


G(0,x, z,T) = P(x,T,Tz), (32) 

rK 1 

AP 0 = e~ rT / G(0, x, T) dz = -e~ rT Q(x, T, TI <), 

Jo 1 

where Q can be obtained numerically by inverting its Laplace transform determined by CD and COD. 

As we shall show in the next section, the integral on the right-hand side of ()31jl decays in the 
order of 0(l/r 2 ), making the numerical integration rather slow. Therefore, we propose a modifica¬ 
tion to improve the efficiency of computation, which turns out to have an economic interpretation 
as well. We introduce 

D(b , x, z, t ) = G(0, x, z, t ) — G(b, x, z, t). 
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It is easy to show that the function 

b —> G(b, x, z, t ) 

is decreasing so D(b, x, z, t ) > 0. Therefore, the price of a conditional Asian option is broke down 
into the price of a regular Asian option and some price spread, 


AP b 


AP 0 + e~ rT 



D{b , x , 2 , T) d z. 


(33) 


Clearly, the regular Asian option price APo can be determined efficiently by methods in Section [2] 
For the computation of D we use 



e st D(b,x, z,t) dt 


1 /'°° 1 

— / —^s^(b,x,s,iTZ,—iT)d,T, 

7T Jo T 


(34) 


where 

4>(6, x , s, a, f3) := F(b, x, s, a, f3) — F{ 0, x, s, a, /3), (35) 

which determines the explicit formula for <!> in (1231) . An advantage of working with D is that the 
integrand of (l34l) decays exponentially as r —>• oo; see Section 6 for details. 


6 Asymptotics as r ^ oo 

The function Y{x) = Y(x,s,iTZ,—ir) satisfies the differential equation 

\g 2 x 2 Y"(x) + rxY'(x) — (s + irz — irx)Y (x) = — 1. (36) 

To get an approximation of Y(x) as r ^ oo, we neglect all terms on the left-hand side of (1361) 
which do not depend on r. We obtain 

1 1 

Y(x, s, irz, — it) ~ —-as r —>• +oo (37) 

i(z — x) T 

and ^ ^ 

— Y(x, s, irz, — ir) ~ —--tt— as r —>• +oo. (38) 

ax i(z — x) z t 

The asymptotics can be proved using |01ver, 1974, Section 10.10]. Numerical tests also confirm 
J3ZD- We can even get the complete asymptotic expansion 

vf ■ • ^ A n (x) 

1 (x,s,trz,-tT) ^ 

n =0 

where A o = as given above and then recursively 

\a 2 x 2 AJ' n + rxA' n — sA n — i(z — x)A n+ 1 = 0. 

Note that (f37l) breaks down when z = x (a turning point appears.) 

With the insight of the asymptotics of Y, we find two numerical issues with the computation 
of the Gurland integral (13TT) . 
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1. Equation (1311) contains the integral 


foo x 

/ —QY(x,s,iTZ,—iT)d,T. (39) 

Jo T 

It follows from (|37l) that the integrand in (1391) behaves like a constant times 1 /r 2 as r —>• oo. 
This decay rate is fast enough to ensure absolute convergence of the integral but the decay is 
very slow which causes problems when one tries to compute the integral (1391) numerically. 

2. When z > x, catastrophic cancellations of the two terms in (1211) may occur in the computation 
of Y(x, s,irz, —it). In general, the computation of Y using (Oil) is numerically stable only 
when z > x. Observe that as r -> oo, each term of the power series of \ iy? in (1211) has a limit, 
namely the series goes term-by-term to the geometric series 



This series converges only when z > x. The result is in agreement with (|37l) . Obviously, the 
second term on the right-hand side of ()21l) must be small compared with the first. However, 
the above is not true when z < x. When we compute Y(b) in (1231) we can assume that z > b 
and there is no computational issue. But when computing the integral (1391) for b < z < x, we 
would face the issue of cancellations. 


Therefore, one should avoid numerical evaluation of (1391) by calculating the price spread in (l33l) 
instead. Now we want to analyze the convergence of the integral in 


foo 1 

/ — 9<3?(6, x, s, irz, — ir) dr. 

Jo T 


(40) 


To determine the behavior of the integrand, we use the following result from fOlve r, 1974, Section 

10 . 8 ] 

/ 7 T \ 1/2 e~ x £ vr 

I Try ) -oTTTI as I ar g A| < —, A —>• oo, 


Kx(Xy) 


\2\J (1 + y 2 ) 1 / 4 


(41) 


where 


C = (l + y 2 ) 1/2 + ln- 


y 


r, % > 0. 


1 + (1 + y 2 ) 1 / 2 ' 

The values y = 0, have to be excluded. We use (fill) with A = 2<r 1 y / 2 Jzt and y\ = —iy/r- If 
z < x we have to use the root with negative imaginary part in the calculation of £i = £(y i). Then 
we obtain 

/ 7T \ !/ 2 e - ^ 1 


F 2 (x, irz,-ir) ~ (-) + ^ )1/4 


as r —y oo. 


Setting y2 = —i\J z anc l C2 = £(2/2) we have 

F 2 (b,s,irz,-i t) ~ (|j) 1/2 (1 ‘J] 1/4 as r -> co. 

From |01ver, 1974. Ex. 7.2, page 378] we obtain in a similar way 

Fi(b) ~ ax 12 o + viYYY 


a 


2A > 


y 2 
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This gives 


F 2 (x) F 2 (x) 

P F 2 (b)-bF'(b)~ bF'(b) 

„ (xy( 1+ M 2 _ V2 _ e -A(a-6) 

\b) v /z 2F6(l + Z/f)i/4(l + y 2 2 )i/4 

Combining this with (137|) . (1381) . we obtain 


3 >( 6 , X, S, ITZ, —ir) r*j 


P 

s 



-(1+m)/2 (7 _ yi _ -A(a-6) 

V-2 irb (1 + 2/i) 1/4 (1 + 2/|) 1/4 


(42) 


It follows from this analysis that x, s, irz, — ir) behaves like a constant times r - 3 / 2 g- a \/r 

as r —> oo with 5Ro > 0. This makes the integral ()40l) numerically more attractive than the one on 
the right-hand side of (13TT) . 


7 The limit a —>■ 0 + 


It is interesting to consider the limit a —> 0 + . We can use this limiting case to gain some insights, 
and also to verify theoretical and numerical results. 

We take a = 0, r < 0, 0 < b < x (if r > 0 then interesting results appear only for x < b.) Then 
the geometric Brownian motion becomes deterministic 

X t = xe rt . 

We define 

T=-Iln|>0. 
r b 

so that xe rT = b. Then 


U t = t AT, 

Vt = — (e r ( tAT ) - 1). 

The joint distribution of (Ut,Vt) is a point mass which moves along a curve in the (u, u)-plane for 
0 < t < T and then stops at t = T. The Laplace transform with respect to (u, v ) is 


[exp i—otUt — &Vt)\ = exp(— a(t A T)) exp(- (e 

r 


_^_( p r(tAT) _ 


m- 


To shorten notation we set 


7 = 


s + a 


Then after some computation we find the triple Laplace transform (HID in the form 


as 


4>o (b,x,s,a,/3)=er( x b) Q 


1 1.7 + h — 

s + a V r 


Y 0 (x,s,a,(3) = —^—M ( 1,7 + 1 ,—) , 
s + a \ r ) 


where M denotes the Kummer function. One can confirm that Yo is the limit of the function (1211) 
as a —> 0 + . In the hyper geometric series 1 F 2 appearing in ()21|) we can go to the limit term by term. 
Obviously, the second term on the right-hand side of (1211) must go to zero. 
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Note that 


x 


z, = h = 


Therefore, 


U t r{t A T) 
G 0 (b,x,z,t ) = 


( e r(iAT) _ ^ 


1 if 7^T) (A<‘ AT » - 1) < z, 

0 otherwise. 


Suppose that 
Then there is a unique t 0 > 0 such that 


-( e ^-1 )<*<*. 


— {e rt0 
rt 0 


1) = z, 


and we have 


s st ° 

e~ st Go(b , x, z, t ) dt =-. 

s 

This allows us to numerically check equation (1311) in the special case a = 0. 
Example: If cr = 0.1, r = —0.2 then 



y(3.0,1.9,2.5.2.1) w 0.094532 
y 0 (3.0,1.9, 2.5, 2.1) « 0.094501 


which is surprisingly close although a is not that small. For one has to make a smaller to get 
good agreement, for example, if a = 0.01, r = —0.2, we get 


<F(2.0,3.0,1.9,2.5.2.1) » 1.873 • 10~ 9 
4>o(2.0,3.0,1.9, 2.5, 2.1) « 1.504 • 10~ 9 


8 Hedging 

Delta-hedging is a common technique used by practitioners to create portfolio for managing in¬ 
vestment risks in over-the-counter derivatives trading. The key element of this dynamic hedging 
strategy is to trade the underlying asset (or futures on the assets) according to the sensitivity 
measure, known as the delta, 

AaP6 = dx ' 

As analytical computation methods of the conditional Asian option is not previously known in the 
literature, BNP Paribas’ calculation of the delta of the conditional Asian option is likely carried out 
by Monte Carlo simulations. Bear in mind that conditional Asian options are long-dated options 
with a typical term of 5 years. A foreseeable difficulty with such an approach is that high efficiency 
of the conditional Asian option price itself can only be achieved at the cost of very expensive 
computations. Since the derivatives are approximated by difference quotients, the computational 
burden is exacerbated with repeated calculations of option prices with small changes in parameters. 
A small sampling error in the estimation of option prices can lead to a large relative error in 
estimating sensitivity measures such as deltas. It may happen that the inevitable sampling errors 
of the estimation statistics overrun subtle differences between option prices resulted from small 
changes in parameters, in which case the estimation of the delta is no longer reliable. Therefore, 
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analytical methods can be of vital importance in the computation of greeks. Although we only 
compute the delta in this paper, the methodology can be easily extended to other greeks. 

Owing to the computational advantage of the price spread, we also break down the delta into 
two parts in the same way as in (f33lh 


A A p 6 = A A p 0 + e rT J ■^~D(b,x,z,T)dz. 
The two quantities can be computed as follows. 


1 


d 


Aap 0 =-e~ rl —Q( x,T,TK), 


where ( d/dx)Q can be determined numerically by its Laplace transform 

P 


9 Q(x,s,y ) = ^ Q 0 


dx 


( 7 ‘ 


4 s cr 2 y 
a 2 ’ Ax 


4s cr 2 y 
a 2 ’ Ax 


Similarly, one can show that ( d/dx)D can be determined numerically by its Laplace transform 


f°° d 1 f°° 1 8 

/ e~ st —D(b,x,z,t)dt = — / — <h(6, x, s, irz, — it) dr, 

Jo ox it J 0 t Ox 


(43) 


where 


d , s p( l -Y(b)) +bY'(b) 

—<h(6, x, s, irz, —it) = ^ F 2 {x)- 


F^(x) = - ^ 1 x~( 3 +^)/ 2 iL A - 


V^ x -(2+ M )/2 / ^ a+i 


a 


a 


-p2(h 


Even though (|43l) works for all z, we can avoid the numerical integration for 0 < z < b by using 

^ D(b,x,z,T ) = -^-G(0,x,z,T) = -^-P(x,T,Tz), 
where ( d/dx)P can be determined numerically by its Laplace transform 


d_ 

dx 


P(x,s,y) = - 


x- 


Po 


4s cr 2 y 
a 2 ’ Ax 


with po given by the explicit expression 


9 Numerical examples 


The computation of regular Asian option price is well-known and extensively studied in the liter¬ 
ature. Nevertheless, we use the first example to confirm the accuracy of the numerical algorithm 
for regular Asian put option price in Section O which serves as part of the numerical algorithm 
for conditional Asian put option price in the second example. Here we use the set of parameters 
under which Asian option prices have been computed in many papers by a variety of techniques 
including Monte Carlo simulations in |Fu et al., 1998], Table 4, Row 2], inverse Laplace transform in 
[Shaw, 2003, Test Problem 2], numerical PDE method in Vecer, 2001] and eigenfunction expansion 
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Fu et al., 1998| 

Shaw, 2003] 

Vecer, 2001 

Linetsky, 2004 

Section [2] 

0.249(MC100) 

0.246416 

0.246 

0.2464156905 

0.2464156819 


Table 1: Verification of regular Asian option prices 


in |Linetsky, 2004, Case #5], etc. Similar comparisons on many more numerical techniques can be 


seen in Fu et_al., 1998] and [Lin etsky , 2004 


r = 0.05, a = 0.50, T= 1, x = 2.0, K = 2.0. 


It should be noted that nearly all works regarding Asian options in the literature are done on call 
options. Here we use the formulas in Section [2] to compute the Asian put option price with the set 
same of parameters first and determine the call option price by the put-call parity. 

In this second example, we consider the pricing and hedging of long-term at-the-money condi¬ 
tional Asian options. Besides demonstrating this new approach of conditional Asian option pricing, 
we also intend to put to the test BNP Paribas’ cost/benefit claim that conditional Asian options 
can achieve 75% of initial delta of regular Asian options with a reduction of 40% in costs. We shall 
use the following set of common parameters 


r = 0.05, b = 1.0, T = 5, x = 2.0, K = 2.0. 


The computation of conditional Asian put option price is carried out in the following five-step 
procedure. 

Step 1 :(D) We calculate the integrals for s = (zln2)/T, i = 1,2, ■ ■ ■ ,2 M, 


r°° i 

/ —$M>(&, x, s, Lrz, —ir) dr (44) 

Jo T 

appearing on the right-hand side of (1M1) for z = b + hk and k = 1,2, ■ ■ ■ , (K — b)/h. The 
integrals are computed by first determining a number a > 0 so that the truncation error of 
(jSP is less than 10 -2 - 2M , 


r i 

D(b,x,z,s)= / —QQ(b,x,s,iTZ,—iT)dT. 

Jo T 

Step 2 :(D) The values of D are computed by the Gaver-Stehfest algorithm 

2 M 

h.D 

n , — _l 

where 


D(b, x, z, T) = 


ln(2) 


k =1 ^ 2 


Hk = (- 1 ) 


M+k 


kf\M 

E 


j M+1 (M\ (2 j\ ( j 


Ml \jj\jj\k-j 


(45) 


i=L(fc+l)/2j 

with l_xj being the integer part of x and k A M = min {k,M}. 

Step 3: (APo) The regular Asian put option price is determined by the inverse Laplace transform method 
outlined in Q, © and (flOl) . 
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Step 4: (Spread) D(b,x,z,T) has been determined in Step 1 for z > b. D(b,x,z,T) = G(0,x, z,T) when 
z < b, which can be determined by (13211 . The quantity Q is then computed numerically by 
inverting its Laplace transform (0 and @ . Then the integral in (1331) is computed last by the 
trapezoidal rule. 

Step 5:(APfr) The final result is produced by (1331) . 


Numerical algorithms in this paper are implemented in Mathematica on a Macbook with 2.9 
GHz Intel Core i5 CPU. The Gaver-Stehfest algorithm for numerical inversion of the Laplace trans¬ 
form utilizes M terms in the Salzer sequence acceleration scheme. It is stated in Abate and Whitt, 2006 


that the algorithm requires a system precision of 2.2 M and produces roughly 0.90M significant dig¬ 
its. In our numerical examples, we use M = 5 and the system precision of 11 digits in Mathematica, 
which is expected to produce about 4 or 5 correct digits for D(b,x, z,t). The integrals (1451) are 
computed with Mathematical numerical integration routine without any specification. In Steps 3 
and 4, G(0, x, z, t ) and f~ G( 0, x, w, t ) dw are computed using the Laplace transform method intro¬ 
duced in Section [2j The needed inverse Laplace transforms are computed with the Euler algorithm 
in A bate a nd Wh itt, 20061- In Step 5, we use the trapezoidal rule with step size of h = 0.1 for the 
numerical integration in (1331) . 

It should be pointed out that only the first step in the five-step procedure, the computation 
of the integrals (1451) , is time consuming whereas all other steps can be done within a second. The 
computation time for the 110 integrals is about 20 minutes. The computation time is longest for z 
between b and x but much shorter when z > x. The reason is that the integrand decays faster as 
r —» oo when z > x. This is confirmed by the asymptotics we derived in Section [6l 

Figure [2] shows the distribution functions of the average-to-date at maturity. The top red line 
represents z —>• P (Yt/T < z) where Yt/T is the (unconditional) average-to-date of asset prices 
with Y defined in Q. The bottom blue line depicts z —>• P (Zt < z) where Z defined in (f28|) 
is the average-to-date above the observation threshold b. The graph provides some geometrical 
interpretation of the option prices. The areas underneath the red line and the blue line between 
z = 0 and z = K represent the accumulated value at maturity of regular Asian put premium and 
conditional Asian put premium (e rT APo and e rT AVf/) respectively. The area below the red line 
and above the blue line accounts for the accumulated value of the price spread between the two 
options. 

We make a comparison of regular and conditional Asian put options in Table [2] under various 
assumptions of volatility coefficient a, which suggests very dramatic changes in the cost/benefit 
comparison of the two types of options. The results are generally not as sensitive to the risk-free 
rate r. While the highlighted case of a = 0.4 does support the claim that the conditional Asian 
put costs 60% of its counterpart Asian put and achieves a delta 75% of that of Asian put, such a 
result should not be over-generalized. This result appears to indicate that the simulations by BNP 
Paribas may have been performed with a volatility at around 40%. 

In general, the difference on per unit cost of delta between a conditional Asian option and 
classical Asian option is much pronounced with high volatility than that with low volatility. For 
example, when a = 0.6, the per unit cost of delta for conditional Asian option is about 0.8675 
whereas that for Asian option is about 1.4389. When a = 0.2, the per unit cost of delta for 
conditional Asian option is about 0.3485 whereas that for Asian option is about 0.3580. An intuitive 
explanation can be given as follows. When the volatility is low, it is less likely that asset prices 
drop below the threshold and hence there is very little difference in average prices between the 
two options. When the volatility is high, it is more likely that asset prices below the threshold are 
filtered out by the average price in conditional Asian option and hence the difference between the 
two options becomes more drastic. 
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Figure 2: Distributions of Z^'- G(0, 2, z, 5) (top) and G(l, 2, z, 5) (bottom). 


a 

ap 6 

AP 0 

APb/APo 

Aap 6 

Aap 0 

Aap^/AaPq 

0.6 

0.1669 

0.4026 

41.47% 

-0.1924 

-0.2798 

68.76% 

0.5 

0.1625 

0.3256 

49.92% 

-0.2029 

-0.2859 

70.97% 








( 0.4 

0.1530 

0.2465 

62.08% 

-0.2156 

-0.2871 

75.11% v . 

0.3 

0.1295 

0.1664 

77.86% 

-0.2316 

-0.2782 

83.27% 

0.2 

0.0810 

0.0877 

92.42% 

-0.2324 

-0.2450 

94.84% 


Table 2: Comparison of regular and conditional Asian options 


Appendix - Lommel’s functions 

The inhomogeneous linear differential equation 

z 2 y" + zy' + (z 2 - A 2 )y = kz^ +1 

has two solutions ks^^\{z) and kS^^z) known as Lommel’s functions, where 

z^+ 1 

= (/r — A + l)(/r + A + 1) 

S^\(z) = s /J , t \(z) 


1-^2 ^1; ^(/x — A + 3), ^(/^ + A + 3); —; 




-l 


: r ~ x + x )) r A + 1 )) [ cos (I(a* - A ) 7r ) J-aW - cos (^(/x + ^) 7r ) A(2)] , 


~^sin(A7r 

and J u is a Bessel function; see [Watson, 1944, Section 10.7]. It is known from |Watson, 1944 
Section 10.75] that 

S/j,,x(z) ~ z^ 1 as z — >• 00 . (46) 

We define a modified Lommel’s function 

T^x(z) = e-^+Ws^iiz) - 2 M_1 T (i(/x - A + 1)) T (i( M + A + 1)) / A (z). 
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Then it can be shown that kT^\(z) is a solution to the differential equation 

z 2 y" + zy' - (z 2 + A 2 )y = kz M+1 . 

The function T is real-valued if both p, A € R and z > 0. Observe that the function can also be 
represented in terms of 


T^\(z) = e e *( a +/x+i)tt /2 cos (i(/x - A)vr) T (i(/x - A + 1)) T (±(/x + A + 1)) ir A (z). 

It follows from (|46l) that 

Tft,\(z) ~ — z^~ l as z —» oo. (47) 

Therefore, the general solution to (11411 is given by 

F(x) = C^x) + C 2 F 2 (x) + , 

where Ci and C 2 are arbitrary constants, and 


k = _2_(VW 

a 2 \ <r 


- 1 - 1 “ 


It follows from the properties of F\,F 2 and (1471) that the solution to (1141) satisfying both boundary 
conditions (1441) and (fT5l) must be F(x) = kx~^ +1 ^ 2 T fJi} x(2^/2/3x / a) which can be rewritten as (12TT) . 
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